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Abstract 

Heavy-tailed distributions are found throughout many naturally occurring 
phenomena. We have reviewed the models of stochastic dynamics that lead 
to heavy-tailed distributions (and power law distributions, in particular) in- 
cluding the multiplicative noise models, the models subjected to the Degree- 
Mass- Action principle (the generalized preferential attachment principle), the 
intermittent behavior occurring in complex physical systems near a bifurcation 
point, queuing systems, and the models of Self-organized criticality. Heavy- 
tailed distributions appear in them as the emergent phenomena sensitive for 
coupling rules essential for the entire dynamics. 

Keywords: Heavy-tailed distributions. Preferential attachment, Intcrmittency, 
Queueing systems, Self-organized criticality 
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1. Introduction 

In 1906, a lecturer in economics at the University of Lausanne in Switzerland 
V.F. Pareto had discovered that the allocation of wealth among individuals can 
be efficiently described by a power law probability distribution, as being self- 
similar over a wide range of wealth magnitudes. A great many such distributions 
have been found in diverse fields of science being a long-known feature of many 
empirical distributions such as Pareto, Zipf, and Levy distributions used to 
model real-world phenomena. In particular, Paul Levy worked on a class of 
probability distributions with " heavy tails" , which he called stable distributions 
largely considered probabilistic curiosities at the time, as heavy-tailed distri- 
butions have properties that are qualitatively different to the many commonly 
used distributions such as exponential, normal or Poisson. 

Since then, there has been a permanent surge of interest to heavy-tailed and 
power law distributions from scientists working in fields as diverse as weather 
forecasting to stock market analysis. The growth rate of the interest to the 
heavy-tailed and power law distributions can be attested by the yearly in- 
crease of the total number of publications on the topic (see Fig. [1]) in the 



arXiv ( http :/ / arxiv. org) , the major forum for disseminating scientific results in 



Physics, Mathematics, Nonlinear Sciences, Computer Science, and Quantitative 
Biology. 

The purpose of this review paper is to explain the models of stochastic 
dynamics that lead to heavy-tailed and power law distributions. 
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Figure 1: The total number of publications devoted to the heavy-tailed and power law distri- 
butions in the arXiv {http://arxiv.org) grows fast each year. 



In this survey, we call a distribution heavy-tailed if it does not have the 
exponentially bounded asymptotes. Namely, given a random variable X char- 
acterized by the probability distribution function 

F{x) = Pt[x<X], 

we say that it has a heavy right tail if 

lim F{x) ■ exp Ax = 00, VA > 0. (1) 

An important subclass of heavy-tailed distributions is the sub- exponential distri- 
butions 01, for which a sum of n independent random variables Xi, X2, . . . , Xn, 
with common distribution F{x) is given. 



Pt[x<Xi+X2 + ...+X„] Pr[x<max(Xi,X2,...,X„)]. (2) 

2:— »oo 

Among sub-exponential distributions, we shall be essentially interested in those 
distributions of a random variable X which are characterized by a power law 
decay, 

Pt[x<X] ^ (3) 

for some 7 > 0. 

There are several physical mechanisms known as underlying the power law 
behaviors. Power law distributions often manifest a form of regularity arising 
through growth processes which is composed of a large number of common 
events and a small number of rarer events happened randomly. As the number 
of events grows, their distribution, under certain conditions, might converge to 
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a steady state ([3]). The mechanism of preferential attachment in which some 
quantity is distributed at random among a number of individuals according to 
how much they already have had been proposed in Q as an algorithm explaining 
power law degree distributions in some scale-free networks. Graphs in that grow 
by successively adding a new vertex say x at each time step that links to an 
already existing one of the degree k, with the probability oc k/N where N is 
the total number of vertices present in the graph. The preferential attachment 
process generates a heavy-tailed distribution following a power law in its tail. 

Power laws are also found in the study of stochastic processes involving 
multiplicative noises [H . A typical equation of multiplicative stochastic process 
is given by a linear Langevin equation with a randomly changing coefficient. 
The effect of such a random coefficient drastically enhance the additive random 
force in the Langevin equation and increase fluctuations. Following let us 
consider temporal fluctuations for a simple discrete time version of the Langevin 
equation, 

x{t + l) - b{t)-xit) + f{t), (4) 

where f{t) represents a random additive noise, and b{t) is a non-negative stochas- 
tic coefficient interpreted as dissipation for b{t) < 1 and amplification for b{t) > 
1. If we assume for simplicity that b{t) and f{t) are independent white noises 
having stationary statistics, and f(t) is symmetric, we obtain, for the second 
order moment, 

{x'{t+l)) = {b') {x\t)) + if) , (5) 

in which the angular brackets denote an average over realizations. For constant 
b^ and there is a stationary solution, 

but (a;^(i)) diverges as i — ?• oo, for (6^) > 1. The statistics of x{t) is estimated 
theoretically by introducing the characteristic function, 

m,t) - (e^«^(*^). (7) 

When (x'^{t)) diverges, t) has singularity at ^ = in the limit oit ^ oo, so 
that Taylor expansion cannot be applied for the steady solution. In such a case 
the following fractional power term can be assumed for the lowest order term 
because the characteristic function is generally a continuous function 

cx)) ^ 1 - const • 1^1'^ + . . . , < /3 < 2, (8) 

which is equivalent to the assumption of power law tails for the cumulative 
probability distribution 

Vr[\x\<X] ^ -j^. (9) 

The validation of stochastic mechanisms generating the power-law behavior re- 
mains an active field of research in many areas of modern science. 



5 



In the forthcoming sections, we consider in details other, more sophisticated 
mechanisms that might bring forth fat-tail statistics in dynamical systems. 

The plan of our review is following. In Sec. [H we discuss the Degree-Mass- 
Action principle in random graphs formation - the preferential attachments and 
their natural generalizations. In Sec. |31 we consider the statistics of bursts 
in systems close to a threshold of instability. Then, in Sec. SI we review the 
emergence of fat tails in queuing systems. Finally, in Sec. [Sj we explain the 
appearance of power law distributions in models of Self-Organized Criticality. 
We conclude in the last section. 



2. Degree-Mass- Action principle in random graphs formation 

Random graphs with scale-free probability degree distributions are ubiqui- 
tous while modeling many real world networks such as the World Wide Web, 
social, linguistic, citation and biochemical networks; an excellent survey is Q. 
The preferential attachment principle, together with its various modifications, 
could be seen as a particular case of degree-mass-action principle since the degree 
of a node acts in that as a positive affinity parameter quantifying the attrac- 
tiveness of the node for new vertices. Our first aim is to construct a family of 
static random graph models, in which vertex degrees are distributed power-law 
like, while edges still have high degree of independence. As usual in random 
graph theory, we will entirely deal with asymptotic properties in the sense that 
the graph size goes to infinity. 

2.1. A random graph space of the preferential attachment model 

We consider graphs with vertex set V = Vn = {l,...,n} where an edge 
between the vertices x and y (denoted by x ^ y) is interpreted as a persistent 
contact between the two nodes. Given x ^V, its degree is denoted by d{x). 

We will think of edges as generated by a pair-formation process in which each 
vertex x - often denoted as an individual - chooses a set of partners according to 
a specified x-dependcnt rule. Therefore the set of individuals which have contact 
with a given vertex x can be divided into two -possible non-disjoint sets: the set 
of nodes which are chosen by x himself and the set of nodes which have chosen 
X as one of their partners. We call the size of the first set the out-degree dout{x) 
of X and the size of the second one the in-degree din{x) of x. Obviously, 

d{x) < dout{x) + din{x), (10) 

and if the choices are sufficiently independent one can expect equality to hold 
almost surely if n — ?■ oo. 

We partition the set of vertices Vn into groups {Ci{ny\^y^ where all members 
of a group Ci{n) choose exactly i partners by themselves {dout = « on Ci{n)). 
Let in,j) the probability for x to choose a fixed partner y G Cj (n) if n 
partners are available for the choice and just one choice will be made be 

Pi(n,j) = A^{n) ^— (11) 
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Here (n) is a normalization constant such that 



A^{n){j2\Q{n)\'-\=l 



and a is a real parameter. Since we want 

WO need \Ci{n)\— to be boimdcd as a fimction of n which will impose later on 
constraints on the constant a. The parameter a acts as an affinity parameter 
tuning the tendency to choose a partner with a high out-degree or low out- 
degrcc. If q: = 0, choices are made without any preferences and Aain) = 1. For 
a > the "highly active" individuals are preferred whereas for a < the "low 
activity" individuals are favored. From this we obtain the basic probability 

Pr[x-y] ^ A^{n)'-^, x e d, yeCj. (12) 

Concerning the size of the sets Ci{n) we will make the following assumption: 

""^^"^1 =P.{n) ^„.oo (13) 



n 

With this choice we have to impose the restriction q < 7 — 1 to ensure the 
convergence of A^ (n). We require furthermore 7 > 2 throughout the paper 
since otherwise the expected in-degree for individuals from a fixed group would 
diverge. The basic probabilities together with the fixed out-degree distribution 
define a probability distribution on each graph with vertex-set Vn, and therefore 
a random graph space Gn (q^iT)- 

First we want to compute the important pairing probabilities. We start with 
the easier case a = 0. 

x-»/ 1 ^ ^ \ i -\~ h ih i -\- k 
PT{x^y\xGCi;yeCk) = , ^n^-cx) (14) 

Likewise one can compute the corresponding probabilities for a 7^ 0. Dropping 
the simple details we just state the result: 

Vy{x y\x (^C^ ]y &Ck) ^ Ao,- n ^ 00. (15) 

It turns out, that the typical graphs in this model still have for a < 2 a power- 
law distribution for the degree with an exponent which can be different from 
the exponent of the out-degree. 

For a > 2 we obtain a degree distribution which follows in mean a power 
law but has gaps. To compare both domains we will use the integrated tail 
distribution 

Fk = Pr(d(a;) > k) . 
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Wc will show that in both cases wc get the same integrated tail distribution. 
Since we are interested in the dependence of the epidemic threshold from the 
power-law exponent of the total degree distribution we have to analyze how this 
exponent varies with the two parameters a and 7 Since the partner choice is 
sufficiently random and not too strongly biased toward high degree individuals 
(that's the meaning of the condition a < 7 — 1) it is easy to see that the in- 
dcgrce distribution of a vertex from group Cj converges for n — > 00 to a Poisson 
distribution with mean const • i". There are essentially two regimes in the 
parameter space, one for which the expected in-degree is of smaller order than 
the out-dcgrce over all groiips and one where the in-dcgrce is asymptotically of 
larger order. In the first case it is clear that the in-degree is too small to have an 
effect on the degree distribution exponent. In other words: the set of individuals 
with degree k consists mainly of individuals whose out-degree is of order k. An 
easy estimation using the formula for the pairing probabilities shows that the 
expected in-degree of individuals from group i is given by 

E {din (x) I x e Ci) ~ const • i" 

asymptotically. Therefore the in-degree is of smaller order than the out-degree 
if a < 1. In the case 7 — 1 > a > 1 the set of individuals with degree k consists 
mainly of individuals from groups with an index of order fc^/". 

2.2. The Cameo Principle. The origins of scale- free graphs in social networks. 

In the present section, we describe an edge formation principle related upon 
a structure apriori imposed on the vertex set. We assume that such a structure 
can be specified by a real positive random variable oj e M+ that quantifies some 
social property of an individual such as its wealth, popularity, or beauty being 
distributed over the population with a given probability distribution (f (w). 

Furthermore, we assume that a link between the two individuals, x and y, 
arises as a result of a directed choice made by either x or y (symbolized by 
X —i- y or y ^ X respectively); in many real life networks edges are formed that 
way. Although the edge creation is certainly a directed process, in the present 
section wc consider the resulting graph to be undirected since for the majority of 
relevant transmission processes defined on the network the original orientation 
of an edge is irrelevant. 

We suppose that the pairing probability follows an inverse mass-action- 
principle: the probability that x decides to connect to y characterized by its 
affinity value u (y) reads as 

Pr{x-^2/|w(y)} ~ aG(0;l) (16) 

where N being the total number of vertices. Let us note that it is not the actual 
value uj (y) which plays a decisive role while pairing, but rather its relative fre- 
quency of appearance over the population. The proposed principle captures the 
essence of antiquity markets - the more rare a property is, the higher is its value, 
and the more attractive it becomes for others. The paring probability model 
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described above had been introduced in [3J] and called the Cameo-principle 
having in mind the attractiveness, rareness and beauty of the small medallion 
with a profiled head in relief called Cameo. And it is exactly their rareness and 
beauty which gives them their high value. 

In the Economics of Location theory introduced by @ and developed by , 
a city, or even more certainly, a particular district in that may specialize in the 
production of a good that can be connected with natural resources, education, 
policy, or just low salary expenditures. City districts compete among themselves 
in a city market not necessarily connected with the quantity of their inhabitants. 
The demand for these products comes into the city district from everywhere 
and can be considered as exogenous. Then, the local attractiveness of a site 
determining the creation of new spaces of motion in that is specified by a real 
positive random variable w > 0. Indeed, it is rather difficult if ever be possible 
to estimate exactly the actual value uj(i) for any site in the urban pattern, since 
such an estimation can be referred to both the economic an cultural factors 
that may vary over the different historical epochs and over the certain groups 
of population. In the framework of a probabilistic approach, it seems natural to 
consider the value as a real positive independent random variable distributed 
over the vertex set of the graph representation of the site uniformly in accordance 
to a smooth monotone decreasing probability density function / (w). 

While investigating the model of Cameo graphs, we assume that 

• The parameter uj is independent identically distributed (i.i.d.) over the 
vertex set with a smooth monotone decreasing density function ip (w) 

• Edges are formed by a sequence of choices. By a choice we mean that 
a vertex x chooses another vertex, say y , to form an edge between y 
and X. A vertex can make several choices. All choices are assumed to be 
independent of each other. 

• If a; makes a choice the probability of choosing y depends only on the 
relative density of ui (y) and is of the form ([T6l) . 

• A pre-defined out-degree distribution determines the number of choices 
made by the vertices. The total number of choices (and therefore the 
number of edges) is assumed to be about const- N. 

We focus on the striking observation that under the above assumptions a 
scale-free degree distribution emerges independently of the particular choice of 
the UJ— distribution. Furthermore, it can be shown that the exponent in the 
degree distribution becomes independent of (p (u) if the tail of if decays faster 
then any power law. 

Let Vn — {1,...,A^} be the vertex set of a random graph space. We are 
mainly interested in the asymptotic properties for N being very large. We 
assign i.i.d. to each element x from the set Vn a continuous positive real random 
variable (r.v.) a; (x) taken from a distribution with density function ip (w). The 
variable uj can be interpreted as a parametrization of Vn ■ For a set 

Cujo^ui = {x : UJ (x) £ [ujo,uJi]} , (17) 
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E(HC^o,-i) = N- U{Lo)duj (18) 



we obtain 



where jlCtjo^t^^ denotes the cardinahty of the set ([T7| . Without loss of generaUty, 
we always assume that > on [0, oo) , and that the the tail of the distribution 
for </j is a monotone function, namely that ip £ C"^ ([0,cxd)) and the second 
derivatives, i'P^) , have no zeros for G (0, /^o) and lj > ujq (/i) . 

Edges are created by a directed process in which the basic events are choices 
made by the vertices. All choices are assumed to be i.i.d. The number of times 
a vertex x makes a choice is itself a random variable which may depend on x. 
We call this r.v. dout (x). The number of times a vertex x was chosen in the 
edge formation process is called the in-degree din (x) . Each choice generates a 
directed edge. We are mostly interested in the corresponding undirected graph. 
If we speak in the following about out-degree and in-degree we refer just to the 
original direction in the edge formation process. Let 

= Fr {x ^ y \ uj ^ uj (y)} 

be the basic probability that a vertex y with a fixed value of a; is chosen by x if 
X is about to make a choice. For a given realization ^ of the r.v. uj over Vjv we 
assume: 

where a G (0, 1) and A (^, N) is a normalization constant. It is easy to see that 
the condition 

oo 





is necessary and sufficient to get 

A{^,N)^A > 0, 

for TV — i> oo where convergence is in the sense of probability. Therefore, we need 
a < 1. 

One might argue that the choice probabilities should depend more explicitly 
on the actual realization ^ of the r.v. uj over Vn -not only via the normalization 
constant. The reason not to do so is twofold. First it is mathematically un- 
pleasant to work with the empirical distribution of lo induced by the realization 
^ since one had to use a somehow artificial N— dependent coarse graining. Sec- 
ond the empirical distribution is not really "observed" by the vertices (having in 
mind for instance individuals in a social network). What seems to be relevant is 
more the common believe about the distribution oi lo. In this sense our setting 
is a natural one. 

The emergence of a power law distribution in the above settings is not a 
surprise as it might seem for the first glance. The situation is best explained by 
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the following example. Let us take 



ip (w) — C ■ e 

and define a new variable 



The new variable w* can be seen as the effective parameter to which the vertex 
choice process applies. What is the induced distribution of ? With 

F{z) = Pr{w* < z} 

we obtain 

F{z) = j ^Hdc. = (20) 



and therefore the w*— distribution 



This is a power law distribution with an exponent depending only on a. 

The detailed results on the degree-degree correlations, the clustering coeffi- 
cients, and the second moment of degree distributions are discussed in j34| . 



3. The statistics of bursts in systems close to a threshold of instability 

Systems driven by random processes at a threshold of stability may exhibit 
a random switching of a signal between a quiescent (stable) and a bursting 
(unstable) state. Such an intermittent behavior is observed over a broad class 
of different systems in physics and nonlinear dynamics. Depending on the origin, 
the intermittent behavior either meets the classification proposed by Pomeau- 
Manneville ^ (the I-III types intermittency) or fits the features of the crisis- 
induced intermittency jlO| . In both cases, the parameters of the models are 
static. Another example of intermittent behavior, called on-ojf intermittency 
has been introduced in [ll| and then observed numerically and experimentally, 
H, E Q E El, [13, EaEi, 0, EH, m ■ The mechanism for this intermittency 



type relies on a random forcing of a bifurcation parameter through a bifurcation 
point. 

The ergodic properties of a system at the threshold of stability can be par- 
tially characterized by the distribution of the quiescent times (the durations of 
laminar phases) P{t) where t G N. Indeed, a complete characterization of the 
statistical properties of the system will imply the knowledge of residence times 
distribution for all the regions of the phase space and not only of the laminar 
regions. However, the former is the first important statistical indicator of such 



11 



dynamics and this is a reason why we focus at the quiescent times distributions 
in the present study. 

Depending on the particular type of intermittency exhibited by the system, 
the statistics of this distribution can asymptoticahy meet either exponential 
laws, or power laws of exponent 7. Particularly, the power-law statistics for 
the quiescent times distribution is claimed to be typical for the systems demon- 
strating the on-off type intermittency, and the value of exponent 7 depends in 
general from the nonlinearity characteristic of the dynamical system considered 
|22| . For example, in the experiments with ion-acoustic instabilities in a labo- 
ratory plasma [23|, due to nonlinear effects, the exponent of power law depends 
on the value of a control parameter. 

In the present section, we discuss the net effect produced on the statistics 
of laminar phases by the stochastic fluctuations of a system state variable (a 
bifurcation parameter) near the fluctuating threshold of stability (a bifurcation 
point). 

We do not refer to any definite physical system displaying an intermittent 
behavior. For the toy model which we introduce, the bifurcation parameter 
and the bifurcation point are considered as random independent variables. It is 
supposed that intermittency takes place in the system when the process crosses 
the threshold value. 

The control parameter of the system is the number 77 G [0, 1], which rep- 
resents a relative frequency of fluctuation of the threshold value: varying the 
parameter 77 amounts to modifying the relation between the characteristic time 
scales of the threshold variable and those of the state variable. 

At = (when the state and the threshold variables have the same time 
scale) the statistics of laminar phases is exponential, while at 77 = 1 (the limiting 
case of quenched threshold) it can be power law; for the intermediate values 
< ry < 1, the statistics is mixed becoming exponential for sufficiently large 
times. 

In general, the statistics of laminar phases depends on the statistics of the 
random system state variable and threshold described by the probability distri- 
butions F and G respectively. For many distributions F and G, the proposed 
model can be solved analytically. 

3.1. Systems at a threshold of instability 

Let us suppose that the state of a system can be characterized by a real 
number x € [0, 1]. Another real number y £ [0, 1] plays the role of a threshold 
of stability. The system is stable as long as x < y and exhibits a sudden 
transition to the irregular state otherwise {x > y). 

We consider Jb 3tS Oil random variable distributed with respect to some given 
probability distribution function 

P{x <u} = F{u). 

In an analogous way, the value of the threshold y is also a random variable 
distributed over the interval [0, 1] with respect to some other probability distri- 
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bution function (pdf) 

P{v<u} = G{u). 

In general, F and G are two arbitrary left-continuous increasing functions sat- 
isfying the normalization conditions 

F(0) = G(0) = 0, = G(l) = 1. 

Given a fixed real number 77 € [0, 1], we define a discrete time random process 
in the following way. At time t = 0, the variable x is chosen with respect to pdf 
F , and y is chosen with respect to pdf G. If x > y, the process continues and 
goes to time t — 1. At time i > 1, the following events happen: 

i) with probability rj, the random variable x is chosen with pdf F but the 
threshold y keeps the value it had at time t — 1. Otherwise, 

ii) with probability 1 — ?7, the random variable x is chosen with pdf _F, and 
the threshold y is chosen with pdf G. 

If a; > J/, the process ends; if a; < y, the process continues and goes to time 
t+l. 

Eventually, at some time step i, when the state variable x exceeds the thresh- 
old value y, the process stops, and the system destabilizes, so this integer value 
t = T acquired in this random process limits the duration of the episode of 
laminar dynamics. In the laps of time, the system regains the composure state, 
when X < y, and the process starts again. 

While studying the above model, we are interested in the distribution of the 
duration of laminar phases Prj{T; F, G) provided the probability distributions F 
and G are given and the control parameter 77 is fixed. 

Even if in our model the state variable x is treated as a random variable, 
what is really important in what follows is the corresponding pdf F. It would 
in fact be possible to treat X 8jS cL deterministic dynamical variable defined by 
the iterated images of a map of the interval [0, 1]. In this case we would assume 
the existence of an invariant ergodic (Bernoulli) measure d-F, for which a; is a 
generic orbit. 

It is also to be noticed that the model proposed resembles closely the coherent- 



noise models [25|, |26| discussed in concern with a standard sandpile model [27 1 
in self-organized criticality, where the statistics of avalanche sizes and durations 
take power-law forms. No exact analytical results concerning the coherent-noise 
models have been obtained so far. The proposed toy model has not been dis- 
cussed in the literature before and, in principle, is much simpler than those 
discussed in 2^ 3] since it does not involve any spatial dynamics typical of 
such extended systems with quenched memory as the original sandpile models. 



3.2. Distribution of residence times below the threshold 

We are interested in the probability Pri{T; F,G) that the random process 
introduced in the previous section ends precisely at time T with a crossing of 
the threshold, provided the distributions F and G are given and rj is fixed. We 
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shall denote Prt{T; F, G) simply as P{T). A straightforward computation shows 
directly from the definitions of Sec. I3.1l that 



P(0) = f dG{y){l~F{v)). (21) 
Jo 

For T > 1, the system can either stay below the threshold in the laminar state 
( a "smvival") {S) or surmount the threshold to a burst state (a "death") [D). 
Both events can take place either in the correlated way (with probability 77; see 
(i) in Sec. 13.11) (we denote them Sc and -De), or in the uncorrelated way (with 
probability 1 — ??; see (ii) in section Sec. 13. ip {Su and £)„). For T = 1, we have 
for example 

P{1) = P[SD,]+P[SD^] 

vJ^dG{y) F{y){l-F{y)) 
+ il-rj)J^dG{y)F{y)J^dG{z)il-F{z)) ^ ' 

7,5(1) + (l-r,)A(l)i?(0). 



Similarly, 

P(2) = 7725(2) + 77(1 -77) (v4(l)S(l) + A(2)B(0)) 

+ (1-77)2^(1)2^(0), 

where where we have defined, for 7i = 0,l,2,..., 



(23) 



A{n) = f dG{y)F{yr (24) 

(25) 



and 

B{n) = dG{y)F{yr{l-F{y)) 
= A{n) - A{7i + 1). 

It is useful to introduce the generating function of P{T): 

00 

P{s) = ^.^P(T). 

The generating property of the function P{s) is such that 

1 d^P(s) 



PiT) 



(26) 



s=0 



Tl dsT^ 

Defining the following auxiliary functions 

p{l) = ri^A{l + l), for l>l, p(0) = 0, 

q{l) ^ (1 -?7)'.4'-l(l), for l>\., q(0) = 0, 

r{l) = [r]B{l + 1) + (1 - ri)A{l + 1)5(0)] , for ^ > 1, r(0) = 0, 
p = r;5(l) + (1-77)^1(1)5(0). 

(27) 
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we find 



Pjs) = B^,)^^,^^ns)+Pmm+pM.msHAil)mr{s)]^ (28) 

1 -p(s)q[s) 

where p(s), g(s), f(s) are the generating functions of p(Z), q(Z), r(/), respectively. 

In the marginal cases of = and t] — I, the probabihty P{T) can be 
readily calculated. For r] = 0, equations (l?fl) and give 

Applying the inverse formula (|26p to equation ((29|). we get 



P,=o(T) = ^^(l)B(O) 



dGiy)F{y) 



T 



dG{v) (1 - F{v)) . 



Therefore, in this case, for any choice of the pdf F{u) and G(m), the probability 
P(T) decays exponentially. 

For ?7 = 1, equations ^ and (gH]) yield 

Pr,=i{s) = B(,s), (30) 

so that ^ 

' dG(y)F(y)^ (1 - F(y)) . (31) 



P,=i(r) = B(T) = 





3.2.1. Some examples of decay in the correlated case rj = 1 

We have just seen that the probability P{T) decays exponentially, in the 
uncorrected case 77 = 0, for any choice of the pdf F and G. 

In the correlated case 77 = 1, many different types of behavior are possible, 
depending on the form of the pdf F and G. We will examine an important class 
of F and G, for which Pri=i{T) can be explicitly computed from equation (|3T|) . 
We will take F and G as absolutely continuous with respect to the Lebesgue 
measure, with 

dF{u) = {1 + a)u°' du, a > — 1, , . 

dG{u) = {1 + (3)il - u)l^ du, /3>-l. ^ 

Here we recognize the family of invariant measures of a map of the interval with 
a fixed neutral point {28l] . 

Equation ([3T|) gives in this case: 

r(2 + /3)r(i + r(i + a)) _ r(2 + /?) r(i + (r + i)(i + a)) 

"-'^ ^ r(2 + /3 + T(l + a)) r(2 + /3 + (T+l)(l + a)) ' 

Using Stirling's approximation, we get for T >> 1: 
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For different values of /3, the exponent of the threshold distribution, we get 
all possible power law decays of P^=i(T). Notice that the exponent (—2 — /3) 
characterizing the decay of Pi,=i(T) is independent of the distribution F of the 
state variable. 

We were not able to prove that the asymptotic decay of Pr^=i{T) is algebraic 
for any choice of the distributions F and G; nevertheless, we have not found 
any counterexample contradicting this conjecture. Let us consider in particular 
the case of uniform F (the results of this section suggest in fact that what 
determines the decay of P{T) is mostly the threshold pdf G): P,^=i{T) is then 
a particular case of a Riemann-Liouville integral, and we did not find any case 
of non-algebraic decay for large T in the tables [29j. 

3.2.2. Upper and lower hounds for P{T) for any rj 
We will use the fact that 

^(1)" < ^(n) < ^(1) and < B{n) < A{1), n = l,2,.... (34) 

The upper bound for A{n) is trivial, since < F{y) < 1 for any y E [0, 1]. The 
lower bound is a consequence of Jensen's inequality, and of the fact that the 
function x — ^ is convex on the interval ]0, 1[ for any integer n. 

We now replace these bounds for A{n) and B{n) in the general formula for 
P(T) and the resulting expressions in all the terms of the sums, except the one 
corresponding to the index n = in Pj(T). (This term, which corresponds to 
a sequence of correlated survivals, has to be treated separately, in order not to 
lose information on the case ry = 1.) Labelling by the index k the number of 
uncorrelated survivals in the sequence of events considered in the sum for P{T), 
we get 

P,{T) < [v^B{T) + ,^^-\l-rj)A{T)B{0)] 

+ Ml) + (1 - v)Ail)Bm YlZl [(1 - ^)^(l)]' V^-'-"" 

(35) 

and 



Pr,{T) > [7fB{T)+if-\\-r,)A{T)B 

+ (1 - v)A[l)B{Q) Y^lZl ll-' [(1 - v)A{irv 



where 7^"^ represents the number of sequences of T — 1 events c, u (c = corre- 
lated survival, u = uncorrelated survival) containing a number k of events u, so 
that 

■ T- 1 
k 



This implies the upper bound 



P,(r) < ifB{T) + {l~r^)A{l)B{Q)[i^+{l-r^)A{l)f 
+ 77^(1) { [77 +(l-r7)A(l)f-i- 77^-1} 



(37) 
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and the lower bound 



Pm > rj^B{T) + {l-^)A{irB{0) 
= 77^P,=i(T) + (l-7y)P,=o(r). 



We thus see that, lor any < 77 < 1, the decay ol distribution P(T) is 
bounded by exponentials. Furthermore, the bounds ((571) a-^d are exact 
in the marginal cases 77 = and r/ = 1. 

3.2.3. Behavior of P{T) for intermediate times 

We have seen in Sec. I3.2.1l that there exists a class ol pdls for which P,-i{T) 
decays like a power law when 77 = 1. In section [3.2.21 we show that for any ry < 1 
the asymptotic decay of PriiT) is exponential. We now make some remarks 
about the behavior of Pr]iT) for 77 close to 1. The first thing to be noted is 
that, for T fixed, Prj{T) is a continuous function of 77, since it is a finite sum 
of continuous functions (see Sec. 13.2.2]) . This results, of course, imply that the 
continuity cannot be uniform in T. This means that, for any fixed interval of 
times [T_ , T+] , with T_ in the range of validity of the power-law asymptotes <\'6'S\i 
of Pi(T), Pti{T) will be arbitrarily close to the same power law for r] sufficiently 
close to 1. For times T T^, the decay becomes exponential. We shall see in 
the next section that for the case of uniform densities, it is possible to estimate 
the value of the crossover time to the exponential behavior as a function of rf. 

3.3. Distribution of quiescent times for the case of uniform densities 

In this section, we consider the distribution of quiescent times for the special 
case of uniform densities dF{u) = dG{u) = du, for all u e [0, 1] and for any 
77 G [0, 1]. In this case, simpler and implicit expressions can be given for P{s) 
and P{T). After some tedious but trivial computation, we get from equation 
HI: 



P{s) = ^ 



l + (l-7;)7(s) 
where 7(5) is defined by 



l + 7(s) 



777(5) 



(39) 



7(.) ^ i^^il^ . (40) 

The asymptotic behavior of P{T) is determined by the singularity of the gen- 
erating function P{s) that is closest to the origin. 

For 77 = 0, the generating function has a simple pole P{s) — [2 ~ s)^^, 
and therefore P{T) decays exponentially, which agrees with the result of the 
previous section. In Fig. 13.31 we have presented the distribution of quiescent 
times P{T) in log-linear scale for 77 = 0. 

For the intermediate values 1 > 77 > 0, the generating function P{s) has 
two singularities. One pole, s — Sq, corresponds to the vanishing denominator 
1 -|- (1 — r])j{s), where sq is the unique nontrivial solution of the equation 

- ln(l - 77s) = s — - — . (41) 
1 - 77 
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2 4 6 8 10 12 14 16 
T 

Figure 2: Distribution of quiescent times for the uniformly distributed variables x and y. 
P-qiT) decays exponentially for r; = 0, consistently with the analytical result P{T) = 2~(^"'"^) 
(solid line). 



Another singularity, s — si — rj~^ , corresponds to the vanishing argument 
of the logarithm. It is easy to see that 1 < sq < si, so that the dominant 
singularity of P{s) is of polar type, and the corresponding decay of P{T) is 
exponential, with rate In (sq (??)), for times much larger than the crossover time 

T,(7;)~ln(so(^))"'. 

The results of Sec. 13.2.21 about the upper bound for the distribution P{T) 
allow us to be more precise about this decay rate. In particular, since B(T) < 
it follows from dST]) that 

P(T) < vMl) + il-v)Ail)BiO) 
which in the case of uniform densities gives 

PiT)<U^Y- (43) 



2 

We have then 

^ (44) 

si So 2 

and we see that the rate In (so(?y)) vanishes like 1 — ?/ as 77 tends to 1. 

When r] tends to 1, the two singularities sq and si merge. More precisely, 
we have 

P,Ms)^ ' + ^'-f^^'--'\ (45) 
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The corresponding dominant term in (pS)) is of order 0(T~^) [30|- This obviously 
agrees with the exact result we get from equation (|3T]) . with dF{u) = dG{u) = 
du : 

P.^.iT) . . (46) 

In Fig. 13.31 we have drawn the distribution of quiescent times P,,=i(T') that 
exhibits the power- law decay, with the slope 7 = — 2. 



-6 



-8 



-10 
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-14 




-16 

3 4 5 6 7 

ln(T) 

Figure 3: Distribution of quiescent times for the uniformly distributed variables x and y. We 
show the power-law decay of P^—i[T) plotted in the log-log scale. The solid line is given by 

In the case of uniform densities, it is possible to get an expression of P{T) 
for all times, and for any value of 77, by applying the inversion formula (|26p to 



P(T) 



(T+l)(T+2) 



'J2k=l {T-k+l)\T~k+2} k J2m=l ( jj'') c{m,k) 

where c(m, k) is defined by 
c(m, k) = ml 



h+t2 + ■■■+l„^=k 
li>l 



(47) 



{h + 1) (fc - ^l) ih + 1) (fc - /l - /a) • • • (Im-l + I)ik-h-l2 Im-l) {U + 1) 
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When 77 ^ 0, there is an alternative way of writing the previous expression: 



^(^) ^ (T+l)(T+2) ^ ^ 

+ J2k=l (T-k+l)(T-k+2) X]i=l(l ~ VY k) , 



where b{l, k) is defined by 



1 

(ii + l)fe + l) ••• {ii-i + l){H + l) ■ 

In Fig. 13. 3[ we have plotted the distribution of quiescent times P,, [T) for the 
intermediate values 77 — 0.5, 77 = 0.7, 77 = 0.9, together with the analytical result 



6(/,A)= ^ 

il+j2H hii = 

ij>0 




Figure 4: Distribution of quiescent times for the uniformly distributed variables x, y at the 
intermediate values r] = 0.5, r] = 0.7, rj = 0.9 (circles). For comparison, the dotted line 2~'^~^ 
presents the exponential decay for = 0, and the dashed line corresponds to 1/(T+ l)(r + 2) 
for T] = 1, see 1146 I I. The solid line is given by I I47II for r] = 0.5, rj = 0.7, rj = 0.9. 



Note that in Fig. 13.31 and Fig. 13.31 (where r/ ^ 1), we only plot distributions 
P(T) up to relatively short quiescent times (T = 16, T = 25), since these times 
are already greater than the crossover time Tc{'q) ^1/ In(so) to the exponential 
decay exp(— In(so)T) (sq defined by equation (|¥T|) ). For much longer times, 
very few survivals are observed, and the statistics gets bad. Of course, Tc(77) 
grows as the parameter 77 tends to 1, so that we have good statistics for longer 
and longer times (in Fig. 13.31 for 77 = 1, the plot is for quiescent times 20 < T < 
2000). 
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A natural question arising in this context is about the relationship between 
the ergodic invariants that quantify the dynamics of deterministic systems, for 
example the Lyapunov exponents, and the scaling laws. The corresponding 
question for models of self-organized criticality is certainly also pertinent since 
in that case a relation is known between the Lyapunov spectrum and the trans- 
port properties [3l|. In our case, however, because of the dynamical character 
not only of the state variable but also of the threshold, some extension of the 
definition of the invariants would be needed, which is beyond the scope of our 
discussion. 



4. Fat tails in queuing systems 



In a simplified model of human activity, 33|, 1471 , |48| , it is viewed as a decision 
based queuing system (QS) where tasks to be executed arrive randomly and 
accumulate before a server S. Under the priority-based scheduling rules, in 
which each incoming task is endowed with a priority index (PI) indicating the 
urgency to process the job, the timing of the tasks follows fat tails probability 
distribution, (i.e the activity of the server exhibits bursts separated by long idle 
periods with the ubiquitous Poisson behavior) . 

There are two types of dynamics: 

i). Service policies based on fixed priority indexes. This case which is consid- 
ered in [S^, I43, IS] assumes that the value a of the PI is fixed once for all. 
Accordingly, very low priority jobs are likely to never be served. To cir- 
cumvent this difficulty [ssj, j43, |48| introduce an ad-hoc probability factor 
< p < 1 in terms of which the limit p 1 corresponds to a deterministic 
scheduling strictly based on the PFs while in the other limit p — >■ the 
purely random scheduling is in use. 

In this setting, the waiting time distribution (WTD) of the tasks before 
service is shown to asymptotically exhibit a fat tail behavior. The main 
point of the Barabasi's contribution is to show that Pl-based scheduling 
rules can alone generate fat tails in the WTD of unprocessed jobs. 



ii). Service policies based on time- dependent priority indexes. Here the prior- 
ity index is time-dependent. This typically models situations where the 
urgency to process a task increases with time and a(t) will hence be rep- 
resented by increasing time functions. 

Clearly, scheduling rules based on such a time-dependent PI do offer new 
specific dynamical features. They are directly relevant in several contexts 
such as: 

a) . Flexible manufacturing systems with limited resource. Here a sin- 
gle server is conceived to process different types of jobs but only a 
single type can be produced at a given time t (i.e. this is the limited 
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resource constraint). Accordingly, the basic problem is to dynam- 
ically schedule the production to optimally match random demand 
arrivals for each types of items. The dynamic scheduling can be opti- 
mally achieved by using time-dependent priority indexes {the Gittins' 
indexes) which specify in real time, which type of production to en- 
gage [40| . Problems of this type belong to the wider class referred as 
the Multi-Armed Bandit Problems in operations research. 

b). Tasks with deadlines. This situation, can be idealized by a queu- 
ing system where each incoming item has a deadline before which it 



definitely must be processed, [43|, |4J, |39|. In this case, to be later 



discussed in the present paper, we can explicitly derive the lead-time 
profile of the waiting jobs obtained under several scheduling rules, 
including the (optimal) time-dependent priority rule known as the 
earliest-deadline-first policy. 

c). Waiting time- dependent feedback queuing systems. In queuing 
networks, priority indexes based on the waiting times can be used 
to schedule the routing through the network. For networks with 
loops, such scheduling policies are able to generate generically stable 
oscillations of the populations contained in the waiting room of the 
queues, [41)]. 

In the context of QS, the waiting time probability distribution (WTD), (i.e. 
the time the tasks spend in the queue before being processed) is a central quan- 
tity to characterize the dynamics. It strongly depends on the arrival and service 
stochastic processes - in particular to the distributions of the inter- arrival and 
service time intervals. The first moments of these distributions, enable to define 
the traffic load 

A 

P = - > 0, 

(i.e. the ratio between the mean service time 1/fi and the mean arrival time 
1/A) and clearly the stability of elementary QS is ensured when < p < 1. 



Focusing on the WTD, [33|, |47|, |48| emphasized that heavy tails in the WTD 



can have several origins, three of which are listed below: 

1) . the heavy traffic load of the server which induces large "bursty" fluctu- 
ations in both the WTD and the busy period (BP) of the QS. For QS with 
feedback control driving the dynamics to heavy trafltic loads, this allows 
to generate self-organized critical (SOC) dynamics, (sj] and the resulting 
fat tails distribution exhibit a decay following a —3/2 exponent. 

2) . the presence of fat tails in the service time distribution produce fat 
tails of the WTD a property which is here independent of the scheduling 
rule [36| . For the convenience of the reader, we give here a short review 
of these results. 

3) . priority index scheduling rules as discussed in 47, 48l |. 
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In this section we pay the essential attention to the case iii) but contrary 
to the discussion carried in H, 47, we shall here consider the dynamics 



in presence of age-dependent priority indexes. As it could have been expected, 
these aging mechanisms generate new behaviors that will be explicitly discussed 
for two classes of models. 

4.1. Waiting time distributions for queuing systems with fat tail service times 
Let us reproduce here a result recently derived in [361] . 

Theorem 4.1 (Boxma). Assume that the (random) service time in a M/G/1 
QS is drawn from a PDF with a regularly varying tail at infinity with index 

V £ (—1, —2), (regularly varying with index v £ —2) =^ fat tail with index 

V G (— 1,— 2)J. For this range of asymptotic behaviors of the PDF, the first 
moment (3 of the service exists. 

Assume further that the service is delivered according to a random order 
discipline. Then the waiting time distribution Wros exhibits a fat tail with 
index (1 — i^) G (—1,0) and more precisely, we can write 

PriWnos >x) ^ C ^''^ ^(^)' (49) 

where p < 1 is the traffic intensity, j3 the average service time, C{x) a slowly 
varying function and 

KPi^) = / fW,P,v)du, 







with: 



The fat tail behavior given in (|49p is therefore entirely inherited from the 
fat tail behavior of the service and is not affected by any reduction of the trafhc 
intensity p. Note also that change of the scheduling rule cannot get rid of this 
fat tail behavior. This point can be explicitly observed in [s^, 5^1, who show 
that for the previous M/G/1 QS with a random order service (ROS) service 
discipline, one obtains: 

Pr(WRos > x) (Xx^oo h{p, v) ■ Pr(VKFCFS > x), (50) 

from which we directly observe that the fat tail in the asymptotic behavior in 
not altered by a change of the scheduling rule. 

Note finally that for the M/M/1 QS, (i.e. exponential service distributions 
and hence no fat tail) , fi^] shows that the random order service scheduling rule 
yields: 

Pr(W^ROS>x) 0C,_,oo CpX-5/6g-7--faV3^ ^^^^ 
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with 



C{p) 



22/33-1/2^5/6^17/12 



[1 




■3 exp 



{ 




and 



7 




P 



1/6 



which has to be compared with the FCFS scheduhng disciphne, which for the 
same M/M/l QS reads as, [3^: 



While the detailed behaviors given by (ISTI) and (15^ clearly differ, they however 
both share, in accord with 133', an exponential decay. 

4-2. Scheduling based on time dependent priority indexes 

The most naive approach to discuss the dynamics of QS with scheduling 
based on time- dependent priority indexes is to think of a population model in 
which the members suffer aging mechanisms which ultimately will kill them. 

Naively, we may consider the population of a city in which members are 
either born in the city or immigrate into it at a certain age and finally die in the 
city. Assuming that the death probability depends on each individual age, the 
study of the age structure of the population exhibits some of the salient features 
of our original QS. 

First, we discuss this class of models and then return to the original model of 
L. Barabasi fss'] to consider a simple QS where each task waiting to be processed 
carries a deadline (playing the role of a PI) and as time flows the these deadlines 
steadily reduce - this implies a {time dependence of the PI). At a given time, 
the scheduling of the tasks follows the " earliest-deadline- first" (EDF) policy and 
given a queue length configuration, we shall discuss the lead-time (lead-time = 
deadline - current time) profile of the tasks waiting to be served. 

4-2.1. Tasks population dynamics with time dependent priority indexes 

Consider a population of tasks waiting to be processed by S with the fol- 
lowing characteristics: 

i) . An infiow of new tasks steadily enters into the queuing system. Each 
tasks is endowed with a priority index (PI) a which indicates its degree 
of urgency to be processed. In general, the tasks are heterogenous as the 
PI are different. In the time interval [t,t + At], the number of incoming 
jobs exhibiting an initial PI in the interval [a, a + Aa] is characterized by 
G(a,i)AtAa. 



Pr(W'FCFS > x) 



1 



(1 - p)e-(i-'')^/^. 



(52) 



24 



ii). Contrary to the situations discussed in [s^, an "aging" process directly 
affects the urgency to process a given task. In other words the priority 
index a is not frozen in time but a = a{t) monotonously increases with 
time t. For an infinitesimal time increase Ai, in the simplest case we shall 
have 

a{t + M) = a{t) + M. 

Here we slightly generalize this and allow inhomogeneous aging rates writ- 
ten as p{a) > meaning that 

a{t + At) = a{t) + p{a)At. 



iii) . The scheduling policy depends on the PI of the tasks in the queue 
and we will focus on the natural policy "process the highest PI first". 

iv) . at time t, a scalar field M{a, t) counts the number of waiting tasks with 
priority index a. Hence M{a, t)Aa is the number with PI p{a) G [a + Aa]. 
Accordingly, the total workload facing the human server S will be given, 
at time t by: 

POO 

L{t) = / M{a,t)da. (53) 

v) . In the time interval [t,t + At], the server S processes tasks with 
an o-dependent rate ii{a)At. Typically /i(a) could be a monotonously 
increasing function of a. As the service rate /i(a) explicitly depend on the 
PI a, it therefore plays an effective role of service discipline. 

The previous elementary hypotheses imply an evolution in the form; 

M{a+p{a)At,t + At)Aa w M{a,t) Aa - fi{a) M{a,t)AaAt 

+G{a,t) AtAa. 

Dividing by AaAt, we end, in the limits Aa and At 0, with the scalar 
linear field equation: 

—M{a,t)+p{a)—M{a,t) + n{a)M{a,t) = G{a,t). (54) 
ot aa 

It is worth to remark that the dynamics given by (IMl) is closely related to the 
famous McKendrick's age structured population dynamics, [sJl- 
Assuming stationarity for the incoming fiow of tasks (i.e. G{a,t) = Gs{a)), the 
linearity of (j54p enables to explicitly write its stationary solution as: 



M(a) = 7r(a) 



G+l -^^dz 
P{z)'^[z) . 



(55) 
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where 

(56) 



with an integration constant C < oo remaining yet to be determined. Assume 
that the PI attached to the incoming jobs do not exceed a hmiting vahie T, 
namely: 

G{a,t) ^l{a <T)G{a,t) ^ G s{a) = l{a < T) G s{a) , (57) 

where I(a < T) is the indicator function. In other words ((57)) indicates that the 
new coming jobs do not exhibit arbitrarily high Pi's. 
This enables to define: 

= r , dz = r ^f'j dz (58) 

Jo p{z)tt[z) Jo p{z)tt{z) 

and (l55t reads as: 

if a < T 



Mia) = 

I 7r(a) 



^ + Jo p(z)7r(z) 

C + V'(T)] if a>T. 



(59) 



The asymptotic behavior of M{a) for a — >■ oo is entirely due to 7r(a), (the square 
bracket terms are bounded by constants) and therefore (1551) and ([55]) imply: 



M{a) « 7r(a) w <J ' " """"ufal 

a when oc k/a 



when 44ocfca'?-\ g > 0, 

(60) 

In view of ([60]) . the following alternatives occur: 

a) . For g < 0, in (jBO]) . the integral M(z) dz does not exist. In this case 
an ever growing population of tasks accumulates in front of the server and 
the queuing process is exploding. 

b) . For q > 0, a stationary regime exists and in this case the constant 
C < oo in ([55]) can be determined by solving: 



Gs{z)dz^ / M{z)fi{z)dz, (61) 
Jo 

which expresses a global balance between the stationary incoming and out 
going flows of tasks. 



c). For q = which implies that 



/i(a) k 
p(a) a 
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([SO]) produces an exponent-k fat tail distribution for M{a) counting the 
number of waiting tasks with PI a in the system. For T < oo and a oo, 
the fat tail of M{a) is populated by long waiting tasks i.e. those having 
spent more than a — T waiting inside the system before being served. In 
the limiting case, for which 

/i(a) = fi = const 

and p{a) — a (i.e. aging directly proportional to time) which leads to q = 
in (pni) . the density M(a) coincides directly with the WTD for a oo. 



This population model shares several features with the Barabasi's model 
namely: 

a). When a stationary regime exists, the function G's(a) which here plays 



the role of the initial PI distribution in [33|, |47|, |48| , does not affect the tail 
behavior given by (1501) . 

b). The scheduling rule here is implicitly governed by the service rate ^(a) 
which itself depend on time as the PI a = a{t) are time-dependent. Note 
that /^(a) directly influences the asymptotic behavior of (j60p . In particular 
for case c), the tail exponent explicitly depends on /i(a). 

Besides the similarities, we now also point out the important differences between 
the present population model and the model discussed in 



a) the service is not restricted to a single task at a given time (i.e. the 
service resource is not limited). Indeed fi{a) describes an average flow 
of service and hence several tasks can be processed simultaneously - (in 
the city population model the service corresponds to death and several 
individual may die simultaneously). 



b) while the fat tail in [33|, |47|, |48| is entirely due to the scheduling rule and 
therefore occurs even for QS far from traffic saturation, this is not so in 
the population model. Indeed in this last case, fat tails are due to heavy 
traffic loads occurring when the flow of incoming tasks nearly saturates 
the server, (this is implied by g = in (j60p ) - for lower loads occurring 
when g > the fat tail in (pOj) disappears. 

4-2.2. Stochastic dynamics. Real-time queuing dynamics 

In this section we will use the results of the real-time queuing theory (RTQS), 



pioneered in [43| , to explore situations where the incoming jobs have a deadline 
- this problem is already suggested in [sl]. Based on [isl H . Isij and (sij . 
first recall the basic hypotheses and the relevant results of RTQS's. Consider a 
general single server QS with arrival and service being described by independent 
renewal processes with average 1 /A respectively 1 //i and finite variances for both 
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renewal processes. Each incoming task arrives with a random hard time relative 
deadline T) drawn from a PDF G{x) with a density g{x): 

Pr{0 < P < x} = G(x), 

with average (P): 



(P) = / {l-G{x))dx= / xg{x)dx. 
Jq Jo 

At a given time t, we define the lead time C to be given by: 

C = V-t, (62) 

Assume now that the lead time C plays the role of a priority index and the 
service is delivered by using the earliest- deadline- first (EDF) rule with preemp- 
tion (i.e. the server always processes the job with the shortest lead time C). 
Preemption implies that whenever an incoming job exhibits a shorter C than 
the one currently in service, this incoming job is processed before, (i.e. pre- 
empts), the currently engaged task which service is postponed. The EDF rule 
directly corresponds to the deterministic policy (i.e. p = 0,7 = in the original 



Barabasi's contribution [33 1. 



At a given time, one can define a probability distribution corresponding to 
the lead time profile (LTP), 

F{x) = Pr{-cx) <C<x}, 

of the jobs waiting in the QS. The LTP specifies the repartition of tasks having 
a given C at time t. Knowing the queuing population Q at a given time, it 
is shown in [SO^ that for heavy traffic regimes, the LTP can, in a first order 
approximation scheme, expressed by a simple analytical form. Indeed, following 
[39|, let us define the frontier J^{Q) > as the unique solution of the equation 



Q 
A 



{l~G{x))dx, a; e [/,oo) C M+. (63) 

HQ) 



Let us also define the frontier J^{Q) as 



when Q < Q*, 
- < ^^</ when Q>Q*, ^''^ 




with Q* being defined by ^{Q) = I where / defined in 

In [3^, it is shown that two alternative regimes can occur: 

a) Jobs served before deadline. When (V) > Q/X and therefore J-{Q) = 
T{Q), the LTP cumulative distribution F{x) takes the form (see Fig. [5]) 



when x < F{Q), 
F{x) = { (65) 

1 - ^ (/r [1 - G(77)] drj) when < T{Q) < x. 
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f F'(x) =f(x) 




Figure 5: Density of the lead time profile f{x) = dF{x)/dx when T(Q) > 0. 



b) Jobs served after deadline. When F{Q) = - Q/X) < 0, the LTP 
cumulative distribution F{x) takes the form (Fig.[6|) 



0, 



F(x) = 



-A(P) 



Q-X{-D)+XV 

1 -Giv)] dv} , 



when X < {V) - ^ < 0, 
when (V) - ^ <x <l, 
when / < X. 



(66) 




Figure 6: Density of the lead time profile f{x) = dF{x)/dx when T(Q) < 0. 

Remark. The ahernative regimes given by ((65|) and ((66)) can be heuristically 
understood by invoking the Little law which connects the average queue length 
(Q) with the average waiting time (W), [s^, 

(Q) ^ X{W), (67) 

a result independent of the scheduling policy. In view of (|63l) and ([67| . one 
obviously suspects that the LTP strongly depends on the sign of the difference 

(p) - M ^ (p> - {w). 

Intuitively, when (W) exceeds (D), it is expected, in the average, that processed 
jobs will be delivered too late and conversely. While the above heuristic argu- 
ments is strictly valid only for the averages, [39|, were able to show that in heavy 
traffic regimes, it also holds also for the LTP given in (|65|) and ([66)) . 

Assuming that the arriving tasks have positive deadlines, the LTP as given 
by (|65|) and (|66l) imply: 
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a). If the left-hand support of the LTP is negative, then a job entering 
into service is aheady late, (case of ((66)) ) see Fig. |6l 



b) . If the left-hand support of the LTP is positive then a jobs enters into 
service with a positive lead time, (case of (|55|) ) see Fig.[SJ Accordingly, it 
is likely that the tasks will be completed before the deadline expired. 

c) . The critical value Q* = {V)/\ for which F{Q) = 0, corresponds to a 
queue length for which customers are likely to become late. Choosing Q 
exactly to Q* , we cannot expect lateness to disappear completely but for 
Q < Q* lateness will be strongly reduced a behavior clearly confirmed by 



d). For deadline distributions G{x) with fat tails, it follows immediately 
from ([65]) and (|66p that the LTP does possess a fat tail. 

4-2.3. "First come first served" (FCFS) scheduling policies 

Choosing the deadline probability density as g{x) = 5{x), (i.e. zero dead- 
line), the EDF scheduling policy directly coincides with the FCFS rule. For this 
case we have Q* = 0, and implies 



Hence the LTP density is given by ([65]) is merely the uniform probability density 
J7[— Q/A, 0], ([— Q/A, 0] being its support). This expresses the fact that in the 
heavy traffic regime p = A//i w 1, the waiting time behaves as Q x (l//i) « 
Q X (1/A) leading to a LPT linearly growing with Q. For general G{x)^ the LTP 
associated with a FCFS scheduling rule will be given by the convolution of the 
deadline distribution G{x) with the uniform distribution [/[— Q/A,0]. Indeed, 
adding the task deadlines with the time spent in the queue, we recover the tasks 
lead-time. Accordingly, in the heavy traffic regime and for a given queue length 
Q, one explicitly knows the LTP's for both the EDF and the FCFS scheduling 
policies thus enabling to explicitly appreciate their respective characteristics. 
In particular, using (f65|) and (|66|) . one can conclude that for a given queue 
length Q, with the FCFS scheduling rule and the associated LTP F{x) being 
the convolution of G{x) with the U [— Q/A, 0], we obtain 



simulation experiments 32|, [3; 



3' 




when Q < 0, 
when Q > 0. 



(68) 







when X < —Q/\, 



(QM) 



[G(e + Q/A)] d^ 



when — Q/A < x < 0, 



I ^ + + 0/A) - G(e)] dO when x>Q, 



(69) 
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where the constant k reads as: 



Q J-(Q/\) 



The latter equation ahows to emphasize the following features: 

i). When the left-hand support of the deadline distribution G{x) is larger 
than Q/\, the left boundary of the support of F{x) is larger than and 
therefore the jobs experience no delay when entering into service. 



ii). If the left-hand support of G{x) is smaller than Q/A, then it may 
happen that the LTP exhibits a negative left-hand support under the 
FCFS policy and a positive left-hand support under the EDF scheduling 
rule. Hence in this last situation, the FCFS policy would deliver tasks 
with lateness while the EDF tasks will be processed in due time. This 
explicitly confirms intuition that EDF is indeed an efficient policy. It has 
been shown that the EDF scheduling rule is optimal for minimizing the 
number of jobs processed after the deadline (46] . 



iii). If G{x) exhibits a fat tail for a; — !> oo so has the LTP and this 
whatever the scheduling rule in use. This can e directly verified from (|69p 
by studying the LTP density f{x) = dF(x)/dx for x oo, we have: 



A 



G{x + ^)-G{x) 
A 



for 



X oo, 



which when G{x) ^ 1 — x and for Q/\ < const takes the form 



for 



X oo. 



(70) 



Hence, the LTP inherits the fat tail property of G{x) and this even when 
using the optimal EDF scheduling rule. 

Below, we focus on a fully explicit illustration involving the Pareto proba- 
bility distribution 







when X < B, 



G{x) 



l-(f)'""'^ whenf>l, Ol, 
which has no moment of order w — 1 or higher. For w > 2, we have 

\uj - I 



(71) 



{V) = 



B. 
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Using (IMl) with I — B, which imphes 



we obtain 



Q* 



B 



BX 


Q 




1 

2 


B- 



XB 



1/(^-2) 



when Q < 
when Q > ^ . 



(72) 



Using and (|66l) . we can find that the LTP distribution reads as 



Q > 



BX 

uj -2 



F{x) = < 



0, 



BX I B^ 



Q < 



XB 

uj~2 



Fix) = 



Q(w-2) V X , 



BX ( B 



Q(w-2) V X 



(I) 



when X < J-{Q), 
when T{Q) < x < B 

when X > B. 
when X < J-{Q), 
when X > J-{Q). 



(73) 



(74) 

The latter equations describe a fat tail with the exponent uj ~ 2. It is worth 
to mention that ((71)) implies that for w > 2 and for 

Q B 

X ^ uj-2' 

the EDF scheduling policy part of the tasks enter into the service before the 
due date expired. Finally note also, that for w < 2, no moments exists for the 
deadline distribution and hence the theory [s^ cannot be applied directly. We 
conjecture that for these regimes no scheduling rule will be able to deliver tasks 
in due time. 

The results obtained for the LTP, enable us to investigate the asymptotic 
properties of the waiting time distribution. Indeed, assume a heavy traffic 
regime with the EDF scheduling policy. Let us also suppose that for a given 
queue length, some jobs are served too late (i.e., the left boundary of the LTP 
is negative). As under the EDF rule, the more urgent jobs are always served 
first, the waiting times of the last jobs in the queue necessarily exceed their 
deadlines. Therefore, when the deadline distribution exhibits a fat tail, so will 
the WTD distribution. Note that while the EDF policy decreases, compared 
with the FCFS rule, the number of jobs served after their deadline, it cannot 
get rid of the fat tail of the WTD, which is due to the fat tail of G{x). This 
result is fundamentally different from the situation that is valid for the frozen 
in time PI models discussed in 3^ 47, 4^, where the fat tail behavior does not 
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depend on G{x) itself. This can be heuristically understood as, in 33l 47. 48| 



the fat tail is mainly due to the low priority jobs, which, as no aging mechanism 



occurs, are likely to never be served. Note that in [33|, |47|, |48| , stable queuing 
models (i.e., those for which the trafHc) [s^ and fat tails of the WTD disappear 
under a FCFS scheduling rule. Indeed without priority scheduling, the WTD 
always follows an exponential asymptotic decaying behavior. In the presence 
of time-dependent PI, all tasks do finally acquire a high priority and this aging 
mechanism precludes the formation of a fat tail solely due to the scheduling 
rule. Accordingly, in the presence of aging PI, the appearance of WTD with fat 
tails is due to G{x). 

4-.3. The importance of adopting performance scheduling policies 

The results for the LTP derived in the preceding section can be directly 
measured on the actual QS. Consider the queue content of a single-stage QS. 
Assume that at a given time, Q is the observed queue content, and at this 
instant take a snapshot of the lead time associated with each waiting item and 
construct the associated LTP (i.e., the histogram of the observed lead times). In 
heavy traffic regimes (i.e., typically 0.95 < p < I leading to stationary average 
queue lengths p/(l — p)) and under the EDF scheduling policy, the LTP will 
approximately be given by (j65l66p . Actual simulation experiments are reported 
in 4^ 44, 33 and [s^, where an excellent agreement between measured data 



and theory is observed. 

From the human activity viewpoint, the explicit expressions of the LTP 
obtained both for the FIFO and EDS policies show clearly that organizing the 
work scheduling is extremely important. As an illustration, consider a situation 
in which the deadline distribution G{x) follows an exponential law: 

G(x) = 1 - e-"" ^ (P) = -. (75) 

a 

For this situation, we compare two different organization policies: 
1. EDF scheduling policy. Introducing (|75l) into (I63p . we obtain 



HQ) = I when Q<A 



, when Q > 

A ^ a 



(76) 



When J-{Q) > (i.e., the upper line in ([76|). it follows from (|65|) that 

{0 when x < F{Q), 

l-^^l^ -hen ^(Q)<a:<0, (77) 

1 - -^e"""" when < x. 

aQ — 

FIFO scheduling policy. 

With G{x) given by ([75]) . the result given in (|69l) reads as 



0, a;<-S, 
F{x) = ^ Q "QI J ' A--^^"' 



(78) 



a; > 0. 



33 



Comparing ([77]) and (1751) . we conclude that in a heavy traffic regime, for a 
given work load Q, the use of EDF enables us to process tasks in due time with 
a high probability while the naive FIFO policy generates large delays. 

Specifically, when Q < A/a, the EDF policy guarantees that most jobs enter 
into service before the deadline (see (177]) ) and will therefore be served before 
deadline, with a high probability. On the contrary, the FIFO policy result given 
in (|78)) (i.e., obtained for a; = in the last line of (|78)) shows that a proportion 
of 1 - (A [1 - e-"^/^] /aQ) jobs enter the service with delays and will therefore 
be late. 

As far as human resources are concerned, this simple model enables us to 
quantify the importance of adopting performance scheduling policies to respond 
to the burn out - generating challenge: deliver more in less time with fewer 
resources. Along the same lines, one of the key rules to avoid burnout is to learn 
to say no to new incoming tasks if the queue length exceeds a threshold. In our 
modeling framework, the critical threshold does depend closely on the level Q* , 
above which lately served tasks (and hence complaints) are unavoidable. 

5. Power law distributions in Self Organized Criticality 

In the last two decades, a meticulous attention has been drawn to the phe- 
nomenon of self-organized criticality (SOC), a property of dynamical systems 
which have a critical point as an attractor, [23, 0, 0, [5l|. A notable fea- 
ture of these models submitted to a power law statistics is that they have no 
characteristic scales, similarly to the scale invariant systems being in a critical 
state. However, unlike systems tackled by the critical phenomena theory the 
critical state in the SOC models seems to be an attractor of the dynamics and 
seems to be achieved without any tuning of control parameters. It is observed in 
slowly-driven non-equilibrium systems with extended degrees of freedom and a 
high level of nonlinearity. The general idea behind SOC models is very appeal- 
ing. Consider for instance Zhangs sandpile model on Z^, where each site has 
an energy variable which evolves in discrete time-steps according to a simple 
"toppling" rule: If a variable exceeds a threshold value, the excess is distributed 
equally among the neighbors. The neighboring sites may thus turn supercritical 
and the process continues until the excess is " thrown overboard" at the system 
boundary. 

What makes this dynamical rule intriguing is that if the toppling is initiated 
from a highly excited state, then the terminal state (i.e., the state where the 
toppling stops) is not the most stable state, but one of many least-stable, stable 
states. Moreover, the latter state is critical in the sense that further insertion of 
a small excess typically leads to further large-scale events. Using the sandpile 
analogy, such events are referred to as avalanches. 

From the very beginning, large theoretical efforts have been made in order 
to understand a true relation between criticality and self-organized criticality 
both exhibiting a power law behavior [s^ [s^. In particular, the use of various 
renormalization group techniques (RG) which proved their exceptional efficiency 
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in justifying scaling properties in the critical phenomena theory 54 |has been 



in the focus of many studies devoted to the SOC phenomena, [5iii,[53,i 



601 161|.. This still deserves a thorough investigation as a potential candidate 
for the "SOC phenomena theory". 

In the critical phenomena theory, the RG method usually helps to establish 
the long time and large scale asymptotic behavior in infinite systems defined 
by the stochastic differential equations with the Gaussian distributed external 
random force 62| that models random boundary conditions. RG is an effective 
method of studying self-similar scaling behavior in such systems. On the con- 
trary, the majority of models exhibiting SOC phenomena are defined on a finite 
piece £ of a discrete lattice by discrete time dynamical rules [13, . More- 
over, in SOC models the energy is usually dissipated at the open boundaries 
of the lattice piece, while, in the majority of critical phenomena, a quenched 
distribution of absorbing defects through the lattice is imposed. 

The primal goal for the SOC phenomena theory is to investigate the scaling 
properties of SOC models, in particularly, to justify the numerically observed re- 
sults 63,(641 on the power law distribution of avalanche sizes in sandpile models 
introduced by Bak, Tang and Wiesenfeldj_ 27[ . In more general formulation, the 
finite size scaling (FSS) hypothesis 6^ 66| is usually assumed in SOC systems. 



P(x,L)=x ^''T^-^-^ , X = {s, a}. 



(79) 



where V) is the probability distribution of occurrence of an avalanche of a 
given size s (the number of sites involved in a relaxation process), area a, and 
time i, L is the size of the lattice piece. If the FSS Ansatz ([79|) is valid, then the 
dynamical exponents Tx and Gx determine the universality class of the model 
27[- 51|. Within the framework of numerous models, the dynamical exponents 
Tx and Ox are found from the different phenomenological approaches, which 
are loosely related to underlying microscopic models, and therefore some doubt 
remains about the universality of representations such as ([7^ . (67 . 68, 69| . 

To identify the scale invariant dynamics, the real-space RG method had been 
applied to the cellular automaton defined on a 2-dimensional lattice[55, 56 1. This 
approach (Dynamically Driven Renormalization Group (DDRG)) deals with the 
critical properties of the system by introducing in the renormalization equations 
a dynamical steady state condition which assumes non-equilibrium stationary 
statistical weights to be used in the calculation (s^]. The fixed points of scal- 
ing transformations define the dynamical exponents whose value are in a good 
agreement with computer simulation data. Nevertheless, it has been shown (55| 
that the fixed points related to these dynamical exponents prescribed by the 
renormalization group are not accessible form the physical domain of parameter 
values. 

An alternative approach referring to the standard critical phenomena the- 
ory is based on the coarse-graining of microscopic evolution rules for the SOC 
models. In a continuous stochastic partial differential equation related 

to the randomly driven models had been proposed although the threshold na- 
ture of the SOC phenomena was not taken into account. A stochastic partial 
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differential equation subjected to a threshold condition and driven by a Gaus- 
sian distributed external random force acting continuously in time has been 
discussed in [s^, III, HI- 

Therein the external random force introduced into the dynamical equation 
simultaneously models: first, the noise risen due to elimination of microscopic 
degrees of freedom; second, unknown (or undefined) boundary conditions; third, 
a mechanism injecting energy into the system which is supposed to act continu- 
ously in time, breaking the time scale separation, and could provoke avalanches 
to overlap. The threshold condition is taken into account by the Heaviside step 
function 6{x). This step function had been regularized as a limit of continuous 
infinitely differentiable functions and then expanded into power series giving 
rise to an infinite series of nonlinearities in the stochastic differential equation 



57| 



Solutions of this nonlinear partial differential equation could be found by 
iterating in the nonlinearities followed by averaging over the distribution of 
the random force. Then, the long time large scale asymptotic behavior of the 
solutions could be established by means of a dynamic RG procedure in the 
spirit of the dynamic RG-approach (H^]. Particularly, the values of dynamical 
exponents could be found in the form of power series in e = 4 — d. However, due 
to an infinite number of nonlinearities risen in the stochastic differential equation 
by the power expansion of step function, the resulting theory calls for an infinite 
number of charges (coupling constants) and, therefore, cannot be analyzed in 



the framework of the standard dynamic RG scheme. In [57[, only the first two 
nonlinear terms have been kept for the RG analysis of the appropriate stochastic 
problem. All higher order terms had been neglected without any estimation of 
their contributions to the long time large scale asymptotic behavior. Let us note 
that the standard power counting analysis of such a model shows convincingly 
that all these terms are of equal importance for the asymptotic behavior and 
all of them have to be taken into consideration on equal footing. 

It is important to emphasize that the correspondence between the models 
of deterministic dynamics and the above stochastic problem is indeed question- 
able and usually lays beyond the studies. The obvious advantage of such a 
coarse graining approach is that it allows to use the modern critical phenom- 
ena techniques of analysis achieving impressive results on the self-similar scaling 
behavior. 

Here, we present the results 5^ Glj on the long time large scale asymptotic 
behavior for the model based on the nonlinear stochastic dynamics equation de- 
rived from the coarse-graining procedure from the discrete rules of deterministic 
dynamics holding all nonlinear terms in check. This task is highly nontrivial and 
of sufficient interest itself stimulating further developments in modern critical 



phenomena theory [54 



5.1. Coarse-graining of microscopic evolution rules for SO C-models 

Recentl y, t wo randomly driven SOC models proposed by Zhang [i^ and by 
Bak et al. [27| (BTW) have been connected to stochastic differential equations 
[st! ]. For the convenience of the reader, we briefly describe the microscopic rules 
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of these SOC models. Both models are defined on a finite piece of d-dimensional 
lattice £ C Z'' in which any site i € C can store some continuously distributed 
variable Ei usually called energy [sij . 

For Zhang's model, the system is perturbed by a random amount of energy 
6E > At a, randomly chosen site i £ £. Once the value Ei exceeds a given 
threshold value Ec, this site becomes active, and transfers all energy to the 
nearest neighbors. As a result, the neighboring sites can be also activated and 
transfer energy to the next neighbors, etc. until it is absorbed at the open 
boundary dC. The avalanche ends when all sites reached a value of energy 
smaller than Ec- The next energy input into the system occurs only when the 
avalanche has stopped. 

For the BTW model, the amount of energy perturbing the system is fixed 
SE — Ec/q where g is a coordination number, and the amount of energy trans- 
ferred to neighbors from an active site is also fixed at Ec- 

For both models, energy is pumped into the system at the small-scale of 
lattice spacing a and then is transferred to a large scale comparable to the size 
of C and actively dissipated at the open boundaries. 

For each i € C, the microscopic evolution rules can be written in the form 
of a stochastic coupled map lattice (SGML), 

E,{t + I) - E,{t) (80) 

= -Yl i(xE,nm{t) + Ec)e {E^,n{t)) - {xE^{t) + Se)^ {E^{t))] + Q{E, t), 
^ mm 

in which Ei{t) is the exceed of energy over the critical value Ec, X = 1 for 
Zhang's model and x = for BTW. The external noise 

CdE, t) = {5E) ■ (5,,„(,) n [1 - ^ (^J W)] (81) 

acts at a slow time scale, and is present when there are no active sites in the 
lattice. Here n(t) is a random vector pointing the site of the lattice piece £ that 
is perturbed with energy (SE) > 0. 

The dynamics of avalanches governed by the SGML (|80|) evolves infinitely 
fast in comparison with the dynamics of energy feeding. It has been pointed 
out [57| that the SGML (|M|) is invariant under spatial translations, rotations, 



and reflections. Furthermore, if x = Oj the equation is invariant under a parity 
transformation of the order parameter, E —E. The SGML (|80II81|) is supplied 
with the absorbing boundary condition EQc{t) = 0. 

The SGML (|80p can be coarse-grained in order to obtain a continuum stochas- 
tic differential equation [H^ for the effective continuum scalar field E{r,t), 

^fe^ = aoA [(xE{r, t) + Ec) 9 {E{r, t))] + /(r, t) (82) 

where ao > is the only dimensional parameter in the model, [a] ~ L^T~^, 
and depends on the lattice spacing a, the unit time step, and the coordination 
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number q. A is the Laplace operator. The noise f{v,t) is a sum of the mul- 
tipUcative external noise depending on the whole lattice state and the internal 
noise that appears due to the elimination of microscopic degrees of freedom. In 
the continuous model, energy is thought to disappear in those regions of lattice 
where /(r, t) <Q and to arrive at the points for which /(r, t) > 0. In a stationary 
state, these processes are obviously balanced, therefore, {f{r,t)) — 0. 

In [s^], the important effect of dissipation at the open boundaries has not 
been taken into account and a quenched distribution of energy absorbing defects 
has been assumed. Time scale separation of dynamics has been also neglected, 
and the noise f{x), x = r, has been understood just as a quenched Gaussian 
process uncorrelated in space and time with a covariance 

{f{x)f{x'))^2TS''{r-r')6{t-t'), (83) 

which is typical for random walks. In the present section, we use a different 
Ansatz for the covariance {f{x)f{x')) which takes the slow-time scale dynamics 
of the stochastic force f{x) into account (see the next section). 

The continuum stochastic partial differential equation requires a regu- 
larization of the step function at zero. Following j57| . we use 

0(i?)=lim ^ ^, (84) 

as a regularizing function where erf(x) = 7r~^/^ exp[— y^] dy is the error 
function and 51 is the regularization parameter. The reason for this choice of 
regularization procedure is that it allows a power expansion with an infinite 
radius of convergence. 

Developing (j84p in powers of E and substituting the series expansion into 
(|82|) . one obtains the strongly nonlinear stochastic partial differential equation 

TI>1 

where the effective coupling constants take different values depending on the 
model: 

A„o= lim e" i;c0("'(O) + — ^^("-i)(0) , neN (86) 

in which 0^"^ (0) is the n— th order derivative of the regularizing function ([84l) at 
zero. The coefficient A„o becomes formally infinite as e ^ oo, however, the series 
in (|85p converges. In the equation ()85p . we have supplied the parameter and 
the coupling constants A„o with index "0" to distinguish their bare values from 
the renormalized analogs which we shall denote in forthcoming sections simply 
as a and A„ consequently. It has been noted [57] that since 0(2ri+2) j-q-j _ ^ 
even coupling constants vanish for the BTW model, whereas they do not for the 
Zhang's one. The set of coupling constants A„o for both models are identical in 
the limit e — > oo. 
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The equation ([55)1 describes the diffusion of energy in Z'^ issued from a source 
defined by f{x), x = r,t. This e quat ion (up to a minor change of notations) 
has been considered in the work [57[ in the whole space and the important 
effect of dissipation at the open boundaries has not been taken into account. 
Alternatively, a small probability of dissipating an amount of energy Ec/ q has 
been assigned for each site when it topples, instead of transferring it to a certain 
neighbor. This procedure expresses the assumption of random boundaries and 
corresponds to a model defined on an infinite lattice with a dissipation for each 
toppling site. We discuss the possible changes to the critical behavior due to 
the presence of regular absorbing boundary in the section 11. 

We study the long time large scale asymptotic behavior in the system gov- 
erned by the stochastic difi^erential equation ([85]) in the whole space supposing 
that the random force f{x) is Gaussian distributed and characterized by the co- 
variance (see the next section) which goes beyond the "white noise" approxima- 



tion studied in[57[. The introduction of the random force f{x) in (|55|) expresses 



the boundary conditions at the random boundaries [57[. 

5.2. Covariance of random forces 

The introduction of the random force into the equation (j85p phenomenologi- 
cally models a consequence of the elimination of microscopic degrees of freedom 
and, at the same time, the injection of energy into the system. We take the 
time scale separation into account supposing that the dynamics of the slow- 
time scale and fast-time scale components of the random force are essentially 
different. Namely, we suppose that in the slow-time scale (i.e., the time scale 
of energy injection) this dynamics can be taken as the white noise like in [st} , 
{F{x)) = 0, with the covariance 

Dp = {F{x)Fix')) = 2T5'^{r - r')S{t - t'), (87) 

where F is the Onsager coefficient. However, in the equation ()85l) defining the 
dynamics of relaxation processes evolving in the fast-time scale j57i | , the random 
force /(a;)introduced into r.h.s. has to be also of fast-time scale. We take it as 
the generalized random walks governed by the linear Langevin equation 

^I^+^f{x)=F{x), x^r,t, (88) 

driven by the slow-time scale "white noise" F{x), where the kernel of the pseudo- 
differential operator has the form 

m{k) = poaofc^"^" (89) 

in the Fourier space. Here, the coupling constant po > and the exponent 
2?7 > are related to the reciprocal correlation time at wave number k, tc{k) = 
k'^^~'^ / PqUq. Dimensional considerations show that the coupling constant po is 
related to the characteristic ultra-violet (UV) ultra-violet momentum scale in 
SOC momentum scale A ~ 1/a by po — A^'' and corresponds to microscopic 
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degrees of freedom expelled from the main equation (1551) as a result of the 



coarse-graining of deterministic dynamical rules 57 1 



The exponent rj corresponds clearly to the anomalous diffusion coefficient 



5l|, |82| z = 2(1 — 7y). Let us note that the scaling form of reciprocal correlation 



time tc{k) has interesting connections with the spectrum of Lyapunov exponent 



for the Zhang model. It has been indeed shown [82|] that the Lyapunov expo- 
nents and the corresponding modes relate to the energy transport in the lattice. 
However, the transport in SOC model is anomalous and the transport modes 
correspond to diffusion modes in a non-flat metric given by the probability that 
a site is active. It is remarkable that the Lyapunov spectrum obeys a simple 
finite scaling form, with an universal exponent r < 1, which is directly related 
to rj by the relation 

d 

The parameter po corresponds to the energy injection rate. 

It is believed in the literature that the SOC regime corresponds to the case 
when the injection rate goes to zero, the dissipation rate goes to zero, such that 
the ratio injection/dissipation goes to zero establishing the time scale separation 



55|,l56|. 



In the framework of critical phenomena theory approach to the problems of 
stochastic dynamics (see [77, 83], for example), the model for the random force 
covariance Dp in ([57)1 is chosen under the following reasons: 

i) for the use of the standard quantum-field RG technique, it is important 
that the function Dp have a power- law asymptote at large k; 

ii) since the covariance in (|87|) is static (oc d{t — t')), the required physical 
dimension [{FF)] — L'^T~^ is put on by a suitable combination of the only 
dimensional parameters in the logarithmic theory (ao and fc); 

iii) the "white noise" assumption (|57|) means that Dp oc Const in the Fourier 
space. To meet this requirement, one introduces a regularization parameter 
e > quantifying the deviation form the logarithmic behavior that is similar to 
the well-know e — 4 — d expansion parameter in the critical phenomena theory 

Hi. 

All above requirements are satisfied by the model 

{F{k, uj)F{~k, uj')) = Dp{k) oc alk^-'^-'^' . (90) 
In this model, 2e is completely unrelated to the space dimension d (in contrast 



to the standard critical phenomena approach [5J, |77| , where usually e = 4 — rf) . 
The logarithmic theory corresponds to the value e = 0. 

Finally, the model for the covariance Dp{k) has to be consistent with the 
form of the linear operator (jSQ]) . Namely, from the equation (|88l) , it follows that 
the covariance Df{uj, k) for the pseudo-random force / introduced in the r.h.s. 
of the main equation (jSSp , 



(/(^)/(^')) = / ^ / -^d^fi^^ k) exp Huj{t - t') + ^k{r - r')] , fc ^ |k|, 

(91) 
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is related to Dp{k) as 



Dfiuj, k) = , (92) 



Then, it is natural that the spectral density ol energy injection, 

W{k) = \ J ^Df{u;,k), (93) 

is independent of the correlation time at given wave number, tc{k). This is true 
if one takes Dpik) oc pQk~^^. Eventually, collecting the latter result with the 
previous Ansatz (|90|) . one arrives at the model 

Dp{k) = po«^fc6-''-2s-2,^ (94) 



Both exponents 2ri and 2e in (|94| are the parameters of the double expansion in 
the r] — e plane around the origin r/ = e = 0, with the additional convention that 
£ = 0{rj). The positive amplitude factor po^"^'' is considered as a dimensionless 
coupling constant (i.e., a formal small parameter of the ordinary perturbation 
theory). For the case of random force uncorrelated in space, Dp{k) oc Const, 
the "real" values of e and f] are taken such that 6 — d = 2{r] + e). The similar 
power-law Ansatz for the correlator of random force has been used to model the 
energy pump into the inertial range of fully developed turbulence 84. 83l|. 



The model (|92|) where the function Dp{k) is defined by (|94|) is then more 
realistic and more reach in behavior than the simple "white noise" assumption 



((83|) discussed in the literature before (for example in the work [57[) since it takes 
into account the finite correlation time of energy field set by interactions at a 
level of microscopic degrees of freedom. It has a formal resemblance with the 
models of random walks in random environment with long-range correlations. 
We note that the similar correlator for random force has been discussed for 
the first time in studies devoted to the anomalous scaling of a passive scalar 
advected by the synthetic compressible flow 



The Ansatz ()92|) that we use contains the previously discussed [57[ model 
([83]) as a special case. Indeed, for the rapid-change limit po — ?• oo, the covariance 
([M)) has the form 

Df{uj,k)^^k^-^-^'+^\ (95) 
Po 

and, for e — 77 = 1 — d/2, one arrives at the model ((83|) uncorrelated in space 
and time with F = ao/2po- 

In the opposite limit of "frozen" configuration of the stochastic force, po — J> 
0, the covariance is static (i.e., independent of time argument {t — t') in t- 
representation) , 

Df{uj, k) Tipoalk'^~'^~^5(Lo). (96) 

The latter case obviously corresponds to an external random force acting con- 
tinuously in time. For e = 4 — d, this random force is uncorrelated in space 
(oc 5{y-y')). 
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5.3. An infinite number of critical regimes in SOC- models 

Quantum field theory formulation for SOC models has been introduced and 



studied in [61| following the general approach developed in 78, 79, 76. |84. 



All correlation functions {E{xi) . . . E{xk)) , x = r, t, and response functions 

{5^ [E{x,)...E{xk)]/5f{x\)...6f{x'J) 

expressing the system response for an external perturbation were renormalized 
by subtracting all ultra-violet superficial divergences from Green's functions. 
An infinite number of renormalization constants has been calculated in the one- 
loop approximation in [6l[ . 

Possible scaling regimes of a renormalizable model are associated with the 
infra-red (IR) stable fixed points of the corresponding differential RG equation. 
The coordinates of fixed points {p, ,A„*} in the infinitely dimensional space 
of coupling constants A„ of the stochastic model (1551) are the solutions of the 
equations 

/3p(/3*,A„,) = /3„(p,, A„*) = 0, n=l,2,...,oo (97) 

where /3p.„— functions, the coefficients in the differential RG equation, are some 
rational functions of the parameters p and A„. Any solution of (|97p is an IR- 
attractive (IR-stable) fixed point of the RG equation if the corresponding Jaco- 
bian matrix 

r _ ^ 
dAfe 

is positively defined (i.e., the real parts of all eigenvalues of the matrix Ji^ are 
positive) for small 77>0, e>?7, 0<p<l, where Aq = p and f3i denotes the 
complete set of the /3— functions of the RG-equation. 

It follows from the explicit expressions for an infinite set of /3— functions 
found in 61] that two coordinates of fixed points, Ai* and A2*, can be chosen 
arbitrary, then all other coordinates and A/t*, k > 2 can be found directly 
from the equations (|97l) . Therefore, the RG differential equation for SOC models 
has a two-dimensional surface of fixed points spanned with Ai, and A2* in the 
infinite dimensional space of coupling constants {p, A„}. 

The complete IR-stability analysis for this manifold of fixed points is a 
formidable task. In the limiting case of "white noise" model ([55]) . taking for- 
mally 

-> 00, 77 = 0, 
it is possible to demonstrate that 

Jik « -(n - l)e5ik < 0, 

where 6ik is the Kronecker symbol, so that there are no IR- stable fixed points 
in SOC for zero correlation time at all wave numbers. The time scale separation 
is mandatory for the existence of a critical regime in SOC. 

In the opposite limit of "frozen" configuration of the stochastic force, the 
fixed points of the RG equation have been shown also IR unstable, since 

Jik = -(n - l)eSik < 0. 
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In a general setting, the matrix Jik in such a case has a block triangular form, 
and therefore its eigenvalues coincide with the diagonal elements which can 
be calculated for any /3„. The stability domains are defined by the roots of 
polynomials in p». For instance, the positivity of the first eigenvalues requires 
that 

and 

- 3 - 5p, 

that is true for |p| < 1. As n grows up, polynomials of any large order can appear 
splitting the stability domain into a number of stable and unstable strips. 

It is important to note that in a multi-charge theory, even if the IR-stable 
fixed points of the RG equation exist, the actual trajectory of the system in the 
multi-dimensional (phase) space of coupling constants (in our case, an infinitely 
dimensional space {p,\k\) starting from the given initial values po, Afeo nray 
not achieve any of them. The trajectory can leave the stability domain (in 
the critical phenomena theory, it is usually interpreted as the first order phase 



< 



transition [5J, |77[) breaking the scaling asymptote. 



5.//. Fat tails in SOC-models 

In the IR-stable critical regimes, the Green functions and the response func- 
tions exhibit scaling behavior characterized by the following "critical dimen- 
sions" : 

A[t] = -AH = -2 + 7„* = 2?7-2, 

A[S] = 277 -e, (98) 

A[£:'] = d + e-2ri. 

We have pointed out before that for the random force uncorrelated in space 
{Dp{k) oc Gonst), the "real" value Sr is taken such that 

3 — d/2 = ?; + Sr- 

In two alternative limiting cases, we have 

Er = 1 - d/2 + r] 

(the "white noise" assumption; the system lacks of IR-attractive fixed points) 
and Sr = 4: — d (the "frozen" configuration of the random force). Substituting 
these values into ([98| , we obtain different critical dimensions for the energy field 
A[£^] and the auxiliary filed A[i?'] listed in the Tab. [1] For instance, for the 
critical dimension of the simplest Green function < EE > (w, k), we obtain 

A [<£;£;>] = 2A[E]-d+ A[t], 

and 

A[<£;^>,t] = 2A[E]-d, 



43 



Table 1: Critical dimensions of the fields E and E' depending on the various models for the 
covariancc Dp 



e AlE] A[E'] 

3-rj-d/2 d/2 + 3(r; - 1) d/2 + 3(1 - rj) 
l-d/2 + r] d/2 + Ti-l d/2 + 1 - -q 

i-d 2rj -4 + d 4-2»? 



for its static anafog, 

< EE >st (fc) = (27r)-^ J duj < EE {uj, k). 

The simplest response function < E'E > evaluates the average size of the relax- 
ation process arisen in the system as a reaction for a point-wise perturbation; 
its Fourier transform is the distribution of avalanche size P(s) observed in nu- 
merical experiments. For < E'E >, in the IR-stable critical regime we obtain 

A[{E'E)] = A[E'] 

+A[E]-d + A[t] = -2 + 27]. 

The squared effective radius 

i?2 = J dxx^ {E{x,t)E' (0,0)) (100) 

of the relaxation process at a moment of time t > started at = at the 
origin x = 0, one can find that it scales as 

^ « R''- (101) 

Indeed, since A[R] ~ —1 (by convention), 



dt 



= -2-A[i], 



from where wc have the result ()10ip . The obtained relation is analogous to 
the well-known Richardson's phenomenological law for the diffusion of passive 
admixtures in the ambient turbulent flows ^80l] . 

We have studied the long time large scale asymptotic behavior for the 
strongly nonlinear stochastic problem which relates both to the Zhang and BTW 
models exhibiting the self-organized critical behavior. The proposed model is 
interesting as itself since it is connected to the problem of nonlinear diffusion of 
the chemically active scalar advection in the turbulent flows [86] . The stochastic 
problem has been considered in the bulk, far from a regular boundary. The en- 
ergy dissipation at the open boundaries has not been taken into account, instead 
a quenched distribution of energy absorbing defects has been assumed. 
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We now make several important comments on the further investigations in 
the framework of RG approach to the stochastic differential equations related 
to the SOC models. 

The first comment is on the possible changes for the critical behavior close to 
the regular open boundary exhibiting the absorbing property. The equation (|85p 
can be considered in a half-space z > where the open boundary coincides with 
the z = plane. Then the effect of boundary would be due to the semi-infinite 
geometry of the system: The absence of sites from one- half space {z < 0) changes 
the energy transfer along the surface. Using the critical phenomena analogy 



one can say that the surface does not become critical simultaneously with the 
bulk, but tends to decouple from the rest of the system. Furthermore, the 
pseudo-random force acting at the boundary is to be always negative to ensure 
the complete dissipation of energy, 

/L=o(r,i) <0. (102) 

This is equivalent to introducing the new field /i^ on the surface that provokes 
a perturbation which can spread inside the system. These two effects are in 
competition: If the coordination number q is large, the toppled amount of energy 
Ec/q dissipated at the open boundary is much smaller than that one transferred 
to neighbors. In this case, the perturbation risen due to hj_ close to the boundary 
cannot propagate into bulk. Otherwise, if energy is rather intensively dissipated 
at the boundary than transferred to the neighboring sites, a critical slope can 
appear. 

Another comment is on possible steps beyond the Gaussian approximation. 
Let us note that the study of composite operators of the type {E')"-{x) would 
manage the corrections for the non-Gaussian distributions of random force. We 
also remember that composite operators are important for the definition of finite 
size scaling corrections to the leading RG predicted asymptotes. 

Scaling, renormalization and statistical conservation laws in the Kraichnan 
model of turbulent advection in the context of the renormalization group im- 
proved perturbation theory have been investigated in [SS] . 



6. Conclusions 

Power laws (heavy-tailed) distributions are found throughout many naturally 
occurring phenomena in physics, and efforts to observe and validate them are 
an active area of scientific research. We have considered a number of stochastic 
dynamical models that might generate power law asymptotic distributions. In 
particular, we have reviewed the stochastic processes involving multiplicative 
noise, Degree-Mass- Action principle (generalized preferential attachment prin- 
ciple), the intermittent behavior occurring in more complex physical systems 
near a bifurcation point, some cases of queuing systems, and the models of 
Self-organized criticality. 

These models might be a ground for many natural complex systems. Heavy- 
tailed distributions appear in them as the emergent phenomena sensitive for 
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coupling rules essential for the entire dynamics. Relationships between the 
rules and the power law statistics are strikingly non-linear, as even a small 
perturbation may cause a large effect, a proportional effect, or even no effect at 
all. 
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